% scripts to red and analyze the waveglider data
%The time stamp reported in the "Time" and "Payload" columns in the CSV 
% files are offset by 7 hours. I have not noticed previously that there 
% is a time stamp in the payload string. The payload time is always in UTC, 
%where as the time stamp is in user  time. So when I export data for 
%2013-Jul-22 00:00 to 23:59, that duration is in Denver time. This is 
%certainly a bug in the Liquid Robotics "Export Data" site. They should 
%not change the "TimeStamp" to user's local time. A record line below for example.

%7/22/13 11:00	LROG16 (SO 2009 4) MAG#1	22.7321	-159.39029	*13.203/18:00:00.5 F:035392.313 S:183 D:-000.6m  L0 0466ms Q:99     
%load and process the WaveGlider data from 2013-05-15 to 2013-06-13

[a, b]  = xlsread('/Users/manojnair/projects/wg/wgrawdata/MAG#2MAY31JULY06.xlsx');


fday_min_mag2 = a(:,1) + datenum(1904,1,1); % Excel date for Mac is based on 1904,1,1

%b(1,:) = []; % first row is header

cPayloadDataStr = cell2mat(b(:,4)); % get the payload data

F_min_mag2 = str2num(cPayloadDataStr(:,22:32));

lat_min_mag2 = a(:,3);
lon_min_mag2 = a(:,4);


[a, b]  = xlsread('/Users/manojnair/projects/wg/wgrawdata/MAG#1MAY15JULY06.xlsx');


fday_min_mag1 = a(:,1) + datenum(1904,1,1); % Excel date for Mac is based on 1904,1,1

%b(1,:) = []; % first row is header

cPayloadDataStr = cell2mat(b(:,5)); % get the payload data

F_min_mag1 = str2num(cPayloadDataStr(:,22:32));

lat_min_mag1 = a(:,3);
lon_min_mag1 = a(:,4);


% Get the depth from the HIGH resolution grid (1 GB !)

[X,Y,Z]=grdread2('/Users/manojnair/projects/wg/bathymetry/himbsyn.bathy.v19.grd');

%get latitude and longitude range of MAG#1 and MAG#2
lat_min = min([min(lat_min_mag1) min(lat_min_mag2)]);
lat_max = max([max(lat_min_mag1) max(lat_min_mag2)]);
lon_min = min([min(lon_min_mag1) min(lon_min_mag2)]);
lon_max = max([max(lon_min_mag1) max(lon_min_mag2)]);


LL = Y >=  lat_min -0.1 & Y <= lat_max + 0.1 ; % maximum limits for the Hawaii region
L = X >= lon_min -0.1 & X <= lon_max + 0.1 ;

Xx = X(L);
Yy = Y(LL);
Zz = Z(LL,L);

clear X Y Z;

depth_mag2 = interp2(Xx,Yy,Zz,lon_min_mag2, lat_min_mag2);
depth_mag1 = interp2(Xx,Yy,Zz,lon_min_mag1, lat_min_mag1);

% verify the interpolation
% scatter(lon_min_mag2, lat_min_mag2, ones([1,length(depth_mag2)]), depth_mag2);
% hold on
% scatter(lon_min_mag1, lat_min_mag1, ones([1,length(depth_mag1)]), depth_mag1);

% Make a location map

imagesc(Xx,Yy,Zz);
set(gca,'FontSize',16);
%colormap('gray');
colorbar;
hold on;
plot(lon_min_mag1, lat_min_mag1,'b.','LineWidth',2);
plot(lon_min_mag2, lat_min_mag2,'r.','LineWidth',2);
legend('MAG1','MAG2')
title('Locations on MAG1 and MAG2 on bathymetry');

%save /Users/manojnair/projects/wg/MAG1MAG2MAY15JULY06

%% Read HON data for May 15 to June 14 in minite files

S = dir('/Users/manojnair/projects/wg/HONdata/*.min');

F_HON = nan([1,length(S)*1440]);
fday_HON = nan([1,length(S)*1440]);

st = 1;


for i = 1:length(S),
    
    fname =  ['/Users/manojnair/projects/wg/HONdata/' S(i).name];
    data_stuct_intermag  = intermagnet(fname,'XYZF');
    
    nlength = length(data_stuct_intermag.DATA(:,4));
    F_HON(st:st+nlength-1) =  data_stuct_intermag.DATA(:,4);
    fday_HON(st:st+nlength-1) = data_stuct_intermag.FDAY;
    
    st = st + nlength;
    
    
end;

save /Users/manojnair/projects/wg/HONMAY15JULY06 F_HON fday_HON;


%% process the data for May 22,23, 2013. This day, WG MAG#1 was off Honolulu in waters
% of about 300 m depth.

load  /Users/manojnair/projects/wg/MAG1MAG2MAY15JULY06 fday_min_mag1 lon_min_mag1 lat_min_mag1 F_min_mag1 depth_mag1
load /Users/manojnair/projects/wg/HONMAY15JULY06.mat  F_HON fday_HON;

% start_date = datenum(2013,5,22);
% end_date = datenum(2013,5,24);

start_date  =  735376.5783410139; % The big region
end_date =    735377.6013824885 ;

start_date = 735376.6569700461;% small region 1
end_date = 735376.8125000000;

% start_date =   735376.8953035275; % small region 2
% end_date =    735377.1126948318;

L = fday_min_mag1 > start_date & fday_min_mag1 < end_date;

fday_HON = fday_HON - 7/24;

LL = fday_HON > start_date & fday_HON < end_date;

lat = lat_min_mag1(L);
lon = lon_min_mag1(L);
F = F_min_mag1(L);
fday = fday_min_mag1(L);


% fit a plane to the data

p = fit([lat lon],F,'poly22') ;

% p = fit([lat lon],F,'poly22','robust','LAR') ;

% Evaluate the model at the observation points

Z = feval(p,lat,lon);
%find the grdient of the plane

d1 = feval(p,0,0);
d2 = feval(p,0,1);
d3 = feval(p,1,0);

dx = (d1-d3)/deg2km(distance(0,0,0,1)); % dF/dx in nT/Km
dy = (d1-d2)/deg2km(distance(0,0,1,0)); % dF/dy in nT/Km

S = sqrt(dx.^2 + dy.^2); % max slope



% now plot all the shit
plot(fday, F-nanmean(F),'r', 'LineWidth',2);
set(gca,'FontSize',16);
hold on;
plot(fday, F-Z,'b', 'LineWidth',2);
plot(fday_HON(LL), F_HON(LL) - nanmean (F_HON(LL)),'g', 'LineWidth',2);
title(datestr(mean(fday)));
datetick;
ylabel('\Delta F nT');
legend('Observed','Corrected','HON');


%% process the data for May 18,19, 2013. This day, WG MAG#1 was off Honolulu in waters
% of about 300 m depth.

load  /Users/manojnair/projects/wg/MAG1MAG2MAY15JULY06 fday_min_mag1 lon_min_mag1 lat_min_mag1 F_min_mag1 depth_mag1
load /Users/manojnair/projects/wg/HONMAY15JULY06 F_HON fday_HON;
L = fday_min_mag1 > datenum(2013,5,18) & fday_min_mag1 < datenum(2013,5,20);

fday_HON = fday_HON - 7/24;

LL = fday_HON > datenum(2013,5,18) & fday_HON < datenum(2013,5,20);

lat = lat_min_mag1(L);
lon = lon_min_mag1(L);
F = F_min_mag1(L);
fday = fday_min_mag1(L);
depth = nanmean(depth_mag1(L));

% fit a plane to the data

p = fit([lat lon],F,'poly22') ;

%p = fit([lat lon],F,'poly44','robust','LAR') ;

% Evaluate the model at the observation points

Z = feval(p,lat,lon);
%find the grdient of the plane

d1 = feval(p,0,0);
d2 = feval(p,0,1);
d3 = feval(p,1,0);

dx = (d1-d3)/deg2km(distance(0,0,0,1)); % dF/dx in nT/Km
dy = (d1-d2)/deg2km(distance(0,0,1,0)); % dF/dy in nT/Km

S = sqrt(dx.^2 + dy.^2); % max slope



% now plot all the shit
plot(fday, F-nanmean(F),'r', 'LineWidth',2);
set(gca,'FontSize',16);
hold on;
plot(fday, F-Z,'b', 'LineWidth',2);
plot(fday_HON(LL), F_HON(LL) - nanmean (F_HON(LL)),'g', 'LineWidth',2);
title(datestr(mean(fday)));
datetick;
ylabel('\Delta F nT');
legend('Observed','Corrected','HON');

%%
% there are some un removed fluctuations still. It could be an aliasing
% issue. Now looking at the seconds data for the same day

% S = dir('/Users/manojnair/projects/wg/seconds/*.xlsx');
%
% lat = nan([1,3600*24]);
% lon = nan([1,3600*24]);
% fday = nan([1,3600*24]);
% F = nan([1,3600*24]);
%
%
% startI = 1;
%
% for i = 1:length(S),
%
%     [a, b]  = xlsread(['/Users/manojnair/projects/wg/seconds/' S(i).name]);
%     b(1,:) = []; % first row is header
%     this_fday = a(:,1) + datenum(1904,1,1) + a(:,2);
%     this_lat = a(:,4);
%     this_lon = a(:,5);
%     cPayloadDataStr = cell2mat(b(:,7)); % get the payload data
%
% this_F = str2num(cPayloadDataStr(:,22:32));
% % subplot(211);
% % plot(this_fday,this_F);
% % datetick;
% % subplot(212);
% % plot(this_lon,this_lat);
% % title(S(i).name);
% % pause;
%
% lat(startI:startI+length(this_lat)-1) = this_lat;
% F(startI:startI+length(this_F)-1) = this_F;
%
% lon(startI:startI+length(this_lon)-1) = this_lon;
% fday(startI:startI+length(this_fday)-1) = this_fday;
%
% startI = startI + length(this_fday);
%
% display(S(i).name);
% end;
%
% L = isnan(lat) | isnan(lon) | isnan(F) | isnan(fday);
%
% lat(L) = [];
% lon(L) = [];
% F(L) = [];
% fday(L) = [];

% load the data
clear

load /Users/manojnair/projects/wg/wgrawdata/seconds/20130523_seconds_data F lat lon fday ;

% Now fit a plane


p = fit([lat' lon'],F','poly11') ; %21 was good

p = fit([lat' lon'],F','poly11','robust','LAR') ;

% Evaluate the model at the observation points

Z = feval(p,lat,lon);


plot(fday, F - nanmean(F),'r', 'LineWidth',1);
hold on;
plot(fday, F-Z,'b', 'LineWidth',1);
plot(data_stuct_intermag.FDAY, data_stuct_intermag.DATA(:,4) - nanmean (data_stuct_intermag.DATA(:,4)),'g', 'LineWidth',2);

%% another way to fit a plane

x = lon';
y = lat';
z = F'; % Equation of the plane containing
k = fday';
% (x,y,z) points is 2*x+5*y+4*z=3

Xcolv = x(:); % Make X a column vector
Ycolv = y(:); % Make Y a column vector
Zcolv = z(:); % Make Z a column vector
%Kconv = k(:);
Const = ones(size(Xcolv)); % Vector of ones for constant term

%Coefficients = [Xcolv Ycolv Kconv Const]\Zcolv; % Find the coefficients

Coefficients = [Xcolv Ycolv Const]\Zcolv; % Find the coefficients


XCoeff = Coefficients(1); % X coefficient
YCoeff = Coefficients(2); % X coefficient
%KCoeff = Coefficients(3);
CCoeff = Coefficients(3); % constant term
%Using the above variables, 
z = XCoeff * x + YCoeff * y + CCoeff;

L=plot3(x,y,z,'r.'); % Plot the original data points
set(L,'Markersize',2*get(L,'Markersize')) % Making the circle markers larger
set(L,'Markerfacecolor','r') % Filling in the markers
hold on
[xx,yy] = meshgrid(min(lon_min):0.00001:max(lon_min), min(lat_min):0.00001:max(lat_min));

zz = XCoeff * xx + YCoeff * yy + CCoeff;
surf(xx,yy,zz) % Plotting the surface
title(sprintf('Plotting plane z=(%f)*x+(%f)*y+(%f)',XCoeff, YCoeff, CCoeff))

Z = XCoeff * lon + YCoeff * lat + KCoeff * fday + CCoeff;

%%
fname =  '/Users/manojnair/projects/wg/hon20130523vmin.min'
data_stuct_intermag  = intermagnet(fname,'XYZF');


plot(fday, F,'b.-')
hold
plot(fday_min + 7/24 , F_min,'r.-')
datetick('x','keeplimits');
ylabel('F nT');
plot(data_stuct_intermag.FDAY, data_stuct_intermag.DATA(:,4),'g')

legend('Seconds data','Minutes data - time shifted by 7 hours','HON data');

%% Read the sentinal data
clear;
close all

sentinal = load('/Users/manojnair/projects/wg/wgrawdata/Sentinel_ May23_formatted.txt');

load /Users/manojnair/projects/wg/MAG1MAG2MAY15JULY06 fday_min_mag1 lon_min_mag1 lat_min_mag1 F_min_mag1 depth_mag1
loaddatastring = '/Users/manojnair/projects/wg/MAG1Location_9data.mat';
load(loaddatastring);
load /Users/manojnair/projects/wg/HONMAY15JULY06 F_HON fday_HON

lat = lat_min_mag1(L)';
lon = lon_min_mag1(L)';
fday = fday_min_mag1(L);
F = F_min_mag1(L)';

LL = fday_HON >= min(fday) + 7/24 & fday_HON <= max(fday) + 7/24;
honf = F_HON(LL);
honfday = fday_HON(LL);
% interpolate HON data at WG sites
honfint = interp1(honfday, honf, fday + 7/24);
honfint(end) = honfint(end-1);


% interpolate HON data at WG sites
honfint = interp1(honfday, honf, fday + 7/24);
honfint(end) = honfint(end-1);
% interpolate HON data at WG sites
honfint = interp1(honfday, honf, fday + 7/24);
honfint(end) = honfint(end-1);
F_resid = ( F - nanmean(F)) - (honfint - nanmean(honfint))';
LL = isnan(F_resid);

% scatter3(lat,lon,F_resid);
% hold on;
p = fit([lat(~LL)' lon(~LL)'],F_resid(~LL)','poly55');
% plot(p);
Z = feval(p,lat,lon);

sentinal(148720:end,:) = [];
sentinal(1:347,:) = [];
fday_sentinel = datenum(2013,5,22) + datenum(0,0,sentinal(:,1),sentinal(:,2),sentinal(:,3),sentinal(:,4));
Y = runmean(sentinal(:,5),60);

% Get HON data for entire Sentinel 

LL = fday_HON >= min(fday_sentinel) + 7/24 & fday_HON <= max(fday_sentinel) + 7/24;
honf = F_HON(LL);
honfday = fday_HON(LL);

plot(fday_sentinel, Y - nanmean(Y), 'b','LineWidth',1);
hold on;
set(gca,'FontSize',16);
plot(fday + 7/24, (F_resid - Z) + (honfint-nanmean(honfint))' ,'r', 'LineWidth',1);
plot(honfday, honf-nanmean(honf),'g', 'LineWidth',1);
title('2013-May-23');
ylabel('\Delta F (nT)');
xlabel('Time Hours')

%% Fit for MAG data on June 1 - 4, 2013
clear
close all

load  /Users/manojnair/projects/wg/MAG1MAG2MAY15JULY06 fday_min_mag1 lon_min_mag1 lat_min_mag1 F_min_mag1 depth_mag1
load /Users/manojnair/projects/wg/HONMAY15JULY06 F_HON fday_HON



L = fday_min_mag1 > datenum(2013,6,1) & fday_min_mag1 < datenum(2013,6,4,0,0,0);

lat = lat_min_mag1(L)';
lon = lon_min_mag1(L)';
F = F_min_mag1(L)';
fday = fday_min_mag1(L);
depth = depth_mag1(L);

LL = fday_HON >= min(fday) + 7/24 & fday_HON <= max(fday) + 7/24;
honf = F_HON(LL);
honfday = fday_HON(LL);
% interpolate HON data at WG sites
honfint = interp1(honfday, honf, fday + 7/24);
honfint(end) = honfint(end-1);


p = fit([lat' lon'],F','poly11') ; %21 was good


% Evaluate the model at the observation points

Z = feval(p,lat,lon);
figure(1)
plot(fday, F - nanmean(F),'r', 'LineWidth',1);
set(gca,'FontSize',16);
hold on;
plot(fday, F-Z,'b', 'LineWidth',1);
plot(fday, honfint - nanmean(honfint),'g','LineWidth',1);
datetick('x','mmm/dd');
ylabel('\Delta F (nT)');
legend('Measured','Corrected','HON')
%find the grdient of the plane

d1 = feval(p,0,0);
d2 = feval(p,0,1);
d3 = feval(p,1,0);

dx = (d1-d3)/deg2km(distance(0,0,0,1)); % dF/dx in nT/Km
dy = (d1-d2)/deg2km(distance(0,0,1,0)); % dF/dy in nT/Km

S = sqrt(dx.^2 + dy.^2); % max slope

% plot plane fit
figure(2)
scatter3(lat,lon,F);
hold
% set(gca,'FontSize',16);
xlabel('Latitude');
ylabel('Longitude');
zlabel('F (nT)');
hold on
plot(p)

%% Find how many loops are there in the data


base_date = datenum(2013,05,15);

for i = 1:31,
    
    L = fday_min_mag1 > base_date & fday_min_mag1 < base_date + 1;
    LL = fday_min_mag2 > base_date & fday_min_mag2 < base_date + 1;
    
    
    if sum(L) > 0,
        
        figure(1);
        plot(lon_min_mag1(L), lat_min_mag1(L),'r-','LineWidth',1);
        title(datestr(base_date));
    end;
    
    
    if sum(LL) > 0,
        
        figure(2);
        plot(lon_min_mag2(LL), lat_min_mag2(LL),'r-','LineWidth',1);
        title(datestr(base_date));
        
    end;
    
    
    pause;
    close all
    base_date = base_date + 1;
end

%% Select data based on time or loction aggregates and save the index
%MAG1
clear
close all
load  /Users/manojnair/projects/wg/MAG1MAG2MAY15JULY06 fday_min_mag1 lon_min_mag1 lat_min_mag1 F_min_mag1 depth_mag1

close all;
mag_azi1 = azimuth(lat_min_mag1(1:end-1),lon_min_mag1(1:end-1),lat_min_mag1(2:end),lon_min_mag1(2:end));
mag_azi1(end+1) = 0;

h=plot(fday_min_mag1,mag_azi1);
set(gca,'FontSize',16);
set(gcf,'Position', [2         958        1675         140]);
% Set tick every two days
tickDays = floor(min(fday_min_mag1)):1:ceil(max(fday_min_mag1));

% Update the tick locations
set(gca, 'XTick', tickDays);

% Update the tick text
datetick('x', 'dd', 'keepticks')
hold on;
flag = 1;
iLocation = 10;

clr = ['r';'b';'g';'k';'y';'c';'r';'b';'g';'k';'y';'c'];

while (flag == 1)
    
    figure(1);
    pos = ginput(2);
    L = fday_min_mag1 > pos(1,1) & fday_min_mag1 < pos(2,1) ;
    
    figure(2);
    my_color = floor(rand(1)*10)+1;
    
    if sum(L) > 0
        
        plot(lon_min_mag1(L), lat_min_mag1(L),clr(my_color),'LineWidth',1);
        
        hold on;
        
    else
        clf(2);
    end;
    figure(1);
    plot(fday_min_mag1(L),mag_azi1(L),clr(my_color));
    set(gca,'FontSize',16);
    
    options.Interpreter = 'tex';
    options.Default = 'Time';
    choice = questdlg('select','Choice','Time','Location','Quit',options);
    
    if strcmp(choice,'Location'),
        
        figure(2);
        pos = ginput(2);
        L = lat_min_mag1 > pos(1,2) & lat_min_mag1 < pos(2,2) & ...
            lon_min_mag1 > pos(1,1) & lon_min_mag1 < pos(2,1) ;
        
        figure(1);
        clf(1);
        plot(fday_min_mag1,mag_azi1);
        set(gca,'FontSize',16);
        % Update the tick locations
        set(gca, 'XTick', tickDays);
        
        % Update the tick text
        datetick('x', 'dd', 'keepticks')
        hold on;
        plot(fday_min_mag1(L),mag_azi1(L),'r');
        
        figure(2);
        clf;
        plot(lon_min_mag1(L), lat_min_mag1(L),'r','LineWidth',1);
        set(gca,'FontSize',16);
        mean_depth = nanmean(depth_mag1(L));
        ndata = sum(L);
        title(sprintf('depth %4.0f number of data %d', mean_depth, ndata));
        
        options1.Interpreter = 'No';
        options1.Default = 'No';
        choice1 = questdlg('select','Save?','Yes','No','Quit',options1);
        
        if strcmp(choice1,'Yes')
            save(['/Users/manojnair/projects/wg/MAG1Location_' sprintf('%d',iLocation) 'data'],'L','mean_depth','ndata');
            set(1, 'PaperPositionMode', 'auto');
            saveas(1,['/Users/manojnair/projects/wg/MAG1Location_' sprintf('%d',iLocation) 'map'],'jpeg');
            set(2, 'PaperPositionMode', 'auto');
            
            saveas(2,['/Users/manojnair/projects/wg/MAG1Location_' sprintf('%d',iLocation) 'time'],'jpeg');
            iLocation = iLocation + 1;
            
            
        end;
        
    elseif strcmp(choice,'Time'),
        
    figure(1);
    pos = ginput(2);
    L = fday_min_mag1 > pos(1,1) & fday_min_mag1 < pos(2,1) ;
    
    figure(2);
    my_color = floor(rand(1)*10)+1;
    
    if sum(L) > 0
        
        plot(lon_min_mag1(L), lat_min_mag1(L),clr(my_color),'LineWidth',1);
        
        hold on;
        
    else
        clf(2);
    end;
        figure(1);
        clf(1);
        plot(fday_min_mag1,mag_azi1);
        set(gca,'FontSize',16);
        % Update the tick locations
        set(gca, 'XTick', tickDays);
        
        % Update the tick text
        datetick('x', 'dd', 'keepticks')
        hold on;
        plot(fday_min_mag1(L),mag_azi1(L),'r');
        
        figure(2);
        clf;
        plot(lon_min_mag1(L), lat_min_mag1(L),'r','LineWidth',1);
        set(gca,'FontSize',16);
        mean_depth = nanmean(depth_mag1(L));
        ndata = sum(L);
        title(sprintf('depth %4.0f number of data %d', mean_depth, ndata));
        
        options1.Interpreter = 'No';
        options1.Default = 'No';
        choice1 = questdlg('select','Save?','Yes','No','Quit',options1);
        
        if strcmp(choice1,'Yes')
            save(['/Users/manojnair/projects/wg/MAG1Location_' sprintf('%d',iLocation) 'data'],'L','mean_depth','ndata');
            set(1, 'PaperPositionMode', 'auto');
            saveas(1,['/Users/manojnair/projects/wg/MAG1Location_' sprintf('%d',iLocation) 'map'],'jpeg');
            set(2, 'PaperPositionMode', 'auto');
            
            saveas(2,['/Users/manojnair/projects/wg/MAG1Location_' sprintf('%d',iLocation) 'time'],'jpeg');
            iLocation = iLocation + 1;
            
            
        end;
        
    elseif strcmp(choice,'Quit')
        break;
        
    end;
    
    
    
    
    
    
    fprintf('%20.10f %20.10f \n',pos(1,1), pos(2,1));
    
end;

%% Select data based on time or location aggregates and save the index
%MAG2

clear

load  /Users/manojnair/projects/wg/MAG1MAG2MAY15JULY06 fday_min_mag2 lon_min_mag2 lat_min_mag2 F_min_mag2 depth_mag2

close all;
mag_azi1 = azimuth(lat_min_mag2(1:end-1),lon_min_mag2(1:end-1),lat_min_mag2(2:end),lon_min_mag2(2:end));
mag_azi1(end+1) = 0;

h=plot(fday_min_mag2,mag_azi1);
set(gca,'FontSize',16);
set(gcf,'Position', [2         958        1675         140]);
% Set tick every two days
tickDays = floor(min(fday_min_mag2)):1:ceil(max(fday_min_mag2));

% Update the tick locations
set(gca, 'XTick', tickDays);

% Update the tick text
datetick('x', 'dd', 'keepticks')
hold on;
flag = 1;
iLocation = 4;

clr = ['r';'b';'g';'k';'y';'c';'r';'b';'g';'k';'y';'c'];

while (flag == 1)
    
    figure(1);
    pos = ginput(2);
    L = fday_min_mag2 > pos(1,1) & fday_min_mag2 < pos(2,1) ;
    
    figure(2);
    my_color = floor(rand(1)*10)+1;
    
    if sum(L) > 0
        
        plot(lon_min_mag2(L), lat_min_mag2(L),clr(my_color),'LineWidth',1);
        
        hold on;
        
    else
        clf(2);
    end;
    figure(1);
    plot(fday_min_mag2(L),mag_azi1(L),clr(my_color));
    set(gca,'FontSize',16);
    
    options.Interpreter = 'tex';
    options.Default = 'Time';
    choice = questdlg('select','Choice','Time','Location','Quit',options);
    
    if strcmp(choice,'Location'),
        
        figure(2);
        pos = ginput(2);
        L = lat_min_mag2 > pos(1,2) & lat_min_mag2 < pos(2,2) & ...
            lon_min_mag2 > pos(1,1) & lon_min_mag2 < pos(2,1) ;
        
        figure(1);
        clf(1);
        plot(fday_min_mag2,mag_azi1);
        set(gca,'FontSize',16);
        % Update the tick locations
        set(gca, 'XTick', tickDays);
        
        % Update the tick text
        datetick('x', 'dd', 'keepticks')
        hold on;
        plot(fday_min_mag2(L),mag_azi1(L),'r');
        
        figure(2);
        clf;
        plot(lon_min_mag2(L), lat_min_mag2(L),'r.');
        set(gca,'FontSize',16);
        mean_depth = nanmean(depth_mag2(L));
        ndata = sum(L);
        title(sprintf('depth %4.0f number of data %d', mean_depth, ndata));
        
        options1.Interpreter = 'No';
        options1.Default = 'No';
        choice1 = questdlg('select','Save?','Yes','No','Quit',options1);
        
        if strcmp(choice1,'Yes')
            save(['/Users/manojnair/projects/wg/MAG2Location_' sprintf('%d',iLocation) 'data'],'L','mean_depth','ndata');
            set(1, 'PaperPositionMode', 'auto');
            saveas(1,['/Users/manojnair/projects/wg/MAG2Location_' sprintf('%d',iLocation) 'map'],'jpeg');
            set(2, 'PaperPositionMode', 'auto');
            
            saveas(2,['/Users/manojnair/projects/wg/MAG2Location_' sprintf('%d',iLocation) 'time'],'jpeg');
            iLocation = iLocation + 1;
            
            
        end;
        
    elseif strcmp(choice,'Quit')
        break;
        
    end;
    
    
    
    
    
    
    fprintf('%20.10f %20.10f \n',pos(1,1), pos(2,1));
    
end;

%% Fit for MAG2 Location # 3 (6/12 to 6/14) depth 3224
clear
close all

fid = fopen('/Users/manojnair/projects/wg/result_table.txt','at');
%load  /Users/manojnair/projects/wg/MAG1MAG2MAY15JULY06 fday_min_mag2 lon_min_mag2 lat_min_mag2 F_min_mag2 depth_mag2
 load /Users/manojnair/projects/wg/MAG1MAG2MAY15JULY06 fday_min_mag1 lon_min_mag1 lat_min_mag1 F_min_mag1 depth_mag1
%loaddatastring = '/Users/manojnair/projects/wg/MAG1Location_5data.mat';
%plane work 
loaddatastring = '/Users/manojnair/projects/wg/MAG1Location_10data.mat';
% loaddatastring = '/Users/manojnair/projects/wg/MAG2Location_5data.mat';
load(loaddatastring);
load /Users/manojnair/projects/wg/HONMAY15JULY06 F_HON fday_HON

% Get the Ap data
[a,b] = xlsread('/Users/manojnair/projects/wg/apdata/kp_may_june_2013.xlsx');
apdata = reshape(a(:,2:end-1)',[1,length(a)*8]);
dummy = datenum(num2str(a(:,1)),'yyyymmdd');
apfday = reshape(repmat(dummy,[1,8])' + repmat([1.5:3:24]'./24,[1,length(a)]),[1,length(a)*8]);

polysurface = 'poly11';
% polysurface = 'poly55';


lat = lat_min_mag1(L)';
lon = lon_min_mag1(L)';
fday = fday_min_mag1(L);
F = F_min_mag1(L)';
% 
% lat = lat_min_mag2(L)';
% lon = lon_min_mag2(L)';
% F = F_min_mag2(L)';
% fday = fday_min_mag2(L);


% remove noise above 10% of median
% LLL = abs( F - nanmedian(F) ) > 300;
% 
% lat(LLL) = [];
% lon(LLL) = [];
% fday(LLL) = [];
% F(LLL) = [];


[fday, ia] = sort(fday);
lat = lat(ia);
lon = lon(ia);
F = F(ia);


LL = fday_HON >= min(fday) + 7/24 & fday_HON <= max(fday) + 7/24;
honf = F_HON(LL);
honfday = fday_HON(LL);
% interpolate HON data at WG sites
honfint = interp1(honfday, honf, fday + 7/24);
honfint(end) = honfint(end-1);
% mean AP
L1 = apfday >= min(fday) + 7/24 & apfday <= max(fday) + 7/24 ;
mean_ap = nanmean(apdata(L1));

%rms of F

L = isnan(F);
[ stats_F] =    rms( F(~L) - median(F(~L))  );
% LP = abs(F-median(F)) > stats_F.robust_s*3;
% F = interp1(fday(~LP),F(~LP),fday)';


%mean lat/long
mean_lat = median(lat);
mean_lon = median(lon);

% distance from HON observatory 
hondist = deg2km(distance(21.32, -157.9996, mean_lat, mean_lon));

% % plot lat long nd find area
plot(lon,lat,'r.');
x = ginput(4);

xdim = deg2km(distance(x(1,1),x(1,2), x(2,1),x(2,2)))*1000;
ydim = deg2km(distance(x(3,1),x(3,2), x(4,1),x(4,2)))*1000;
close all;

p = fit([lat' lon'],F',polysurface) ; %for MAG1 Location 1

% p = fit([lat' lon'],F','poly11','robust','LAR') ;

% Evaluate the model at the observation points

Z = feval(p,lat,lon);


plot(fday , F - nanmean(F),'r', 'LineWidth',1);
hold on;
plot(fday , F-Z,'b', 'LineWidth',1);
plot(fday, honfint-nanmean(honfint),'g', 'LineWidth',1);
set(gca,'FontSize',16);
datetick('x','mmm/dd');
ylabel('\Delta F (nT)');
legend('Measured','Corrected','HON');
%find the grdient of the plane

d1 = feval(p,0,0);
d2 = feval(p,0,1);
d3 = feval(p,1,0);

dx = (d1-d3)/deg2km(distance(0,0,0,1)); % dF/dx in nT/Km
dy = (d1-d2)/deg2km(distance(0,0,1,0)); % dF/dy in nT/Km

S = sqrt(dx.^2 + dy.^2); % max slope

% find rms

% rms of F

L = isnan(honfint);
[stats_H] =    rms( honfint(~L) - mean(honfint(~L))  );

L = isnan(honfint) | isnan(F)';
[stats_D] =  rms( (honfint(~L) - mean(honfint(~L)) ) - ( F(~L) - mean(F(~L)))'  );

L = isnan(honfint) | isnan(F)';
[stats_Z] =    rms( (honfint(~L) - mean(honfint(~L)) ) - ( F(~L) - Z(~L))'  );

title([loaddatastring([end-21:end-17 end-8]) polysurface ',r1# ' sprintf('%4.2f',stats_D) ',r2# ' sprintf('%4.2f',stats_Z) ...
    ',' sprintf('%4.0fm',abs(mean_depth)) ',d1#m' sprintf('%4.0f',xdim) ',d2#m ' sprintf('%4.0f,Dist %3.0f',ydim, hondist)]);

% duration of meansurements

measurement_duration = max(fday)-min(fday);

% now print the value

fprintf(fid, '%s %s %s %s %10.7f %10.7f %10.2f %10.5f %10.2f %10.7f %10.7f %f %10.5f %10.5f %10.5f %10.5f\n', ...
    loaddatastring([end-21:end-17 end-8]),polysurface, datestr(fday(1)), datestr(fday(end)), mean_lat, mean_lon, ... 
    mean_depth, measurement_duration, S, xdim, ydim, ndata, ...
    stats_F,stats_H,stats_D,stats_Z);

fprintf( '%10.7f %10.7f %10.2f %10.5f %10.2f %10.7f \n %10.7f %d %10.5f %10.5f %10.5f %10.5f\n', ...
     mean_lat, mean_lon, mean_depth, measurement_duration, S, xdim, ydim, ndata, ...
    stats_F,stats_H,stats_D,stats_Z);

% save the figure

set(1, 'PaperPositionMode', 'auto');
saveas(1,['/Users/manojnair/projects/wg/figures/' loaddatastring(end-21:end-8) '_TS'],'jpeg');
% 
% figure(2);
% for i = 1:length(lat),
% pp(i) = median( deg2km(distance(lat(i),lon(i),lat,lon))*1000);
% end;
% hist(pp);
% xlabel('Dist (m)');
% ylabel('ndata');
% title([loaddatastring([end-21:end-17 end-8]) ',r1# ' sprintf('%4.2f',stats_D) ',r2# ' sprintf('%4.2f',stats_Z) ...
%     ',' sprintf('%4.0fm',abs(mean_depth)) ',d1#m' sprintf('%4.0f',xdim) ',d2# ' sprintf('%4.0fm',ydim)]);
% 
% set(2, 'PaperPositionMode', 'auto');
% saveas(2,['/Users/manojnair/projects/wg/figures/' loaddatastring(end-21:end-8) '_hist'],'jpeg');


%% Select data when MAG1 and MAG2 were together
clear
load  /Users/manojnair/projects/wg/MAG1MAG2MAY15JULY06
load /Users/manojnair/projects/wg/HONMAY15JULY06 F_HON fday_HON

[fday_12,ia,ib] = intersect(fday_min_mag1, fday_min_mag2);

fmag1 = F_min_mag1(ia);
fmag2 = F_min_mag2(ib);
lat1 = lat_min_mag1(ia);
lon1 = lon_min_mag1(ia);
lat2 = lat_min_mag2(ib);
lon2 = lon_min_mag2(ib);
depth1 = depth_mag1(ia);
depth2 = depth_mag2(ib);

LL = fday_HON >= min(fday_12) + 7/24 & fday_HON <= max(fday_12) + 7/24;
honf = F_HON(LL);
honfday = fday_HON(LL);
% interpolate HON data at WG sites
honfint = interp1(honfday, honf, fday_12 + 7/24);
honfint(end) = honfint(end-1);


d = distance(lat1, lon1, lat2, lon2);
mag_azi1 = azimuth(lat1(1:end-1),lon1(1:end-1),lat1(2:end),lon1(2:end));
mag_azi1(end+1) = 0;

mag_azi2 = azimuth(lat2(1:end-1),lon2(1:end-1),lat2(2:end),lon2(2:end));
mag_azi2(end+1) = 0;

h=plot(fday_min_mag1,mag_azi1);


LL = deg2km(d)*1000 < 15000;

plot(fday_12(LL), fmag1(LL));
hold on;
plot(fday_12(LL), fmag2(LL),'r');

rms1 = rms(fmag1(87:419)- fmag2(87:419) ); % manually selected
rms2 = rms(fmag1(35706:36488)- fmag2(35706:36488));

subplot(121);
plot(fday_12(87:419), fmag1(87:419));
set(gca,'FontSize',16);

hold on;
plot(fday_12(87:419), fmag2(87:419),'r');
plot(fday_12(87:419), honfint(87:419),'g');
legend('MAG1','MAG2','HON');
title(sprintf('%s RMS %4.2f',datestr(round(mean(fday_12(87:419)))), rms1))
xlabel('Hours PST');
ylabel(' F (nT)');

datetick;
subplot(122);

plot(fday_12(35706:36488), fmag1(35706:36488),'b');
set(gca,'FontSize',16);
hold on;
plot(fday_12(35706:36488), fmag2(35706:36488),'r');
plot(fday_12(35706:36488), honfint(35706:36488),'g');
legend('MAG1','MAG2','HON');
title(sprintf('%s RMS %4.2f',datestr(round(mean(fday_12(35706:36488)))), rms2))
xlabel('Hours PST');
ylabel(' F (nT)');

datetick;




%% To show the effect of crustal anomaly when the depth is very low
clear;

load /Users/manojnair/projects/wg/MAG1MAG2MAY15JULY06 fday_min_mag1 lon_min_mag1 lat_min_mag1 F_min_mag1 depth_mag1
loaddatastring = '/Users/manojnair/projects/wg/MAG1Location_2data.mat';
load(loaddatastring);
load /Users/manojnair/projects/wg/HONMAY15JULY06 F_HON fday_HON

lat = lat_min_mag1(L)';
lon = lon_min_mag1(L)';
fday = fday_min_mag1(L);
F = F_min_mag1(L)';

LL = fday_HON >= min(fday) + 7/24 & fday_HON <= max(fday) + 7/24;
honf = F_HON(LL);
honfday = fday_HON(LL);
% interpolate HON data at WG sites
honfint = interp1(honfday, honf, fday + 7/24);
honfint(end) = honfint(end-1);


% interpolate HON data at WG sites
honfint = interp1(honfday, honf, fday + 7/24);
honfint(end) = honfint(end-1);
% interpolate HON data at WG sites
honfint = interp1(honfday, honf, fday + 7/24);
honfint(end) = honfint(end-1);
F_resid = ( F - nanmean(F)) - (honfint - nanmean(honfint))';
LL = isnan(F_resid);

scatter3(lat,lon,F_resid);
hold on;
p = fit([lat(~LL)' lon(~LL)'],F_resid(~LL)','poly55');
plot(p);

Z = feval(p,lat,lon);
plot(fday , F_resid,'b', 'LineWidth',1);
hold;
plot(fday , F_resid - Z,'r', 'LineWidth',1);
set(gca,'FontSize',16);
datetick('x','keeplimits');


plot(fday , F - nanmean(F),'r', 'LineWidth',1);
hold on;
plot(fday , F-Z - nanmean(F-Z),'b', 'LineWidth',1);
plot(fday, honfint - nanmean(honfint),'g', 'LineWidth',1);
set(gca,'FontSize',16);
datetick('x','mmm/dd');
ylabel('\Delta F (nT)');
legend('Measured','Corrected','HON');


%% May 23 curve processin. To show the effect of crustal anomaly when the depth is very low
clear;

load /Users/manojnair/projects/wg/MAG1MAG2MAY15JULY06 fday_min_mag1 lon_min_mag1 lat_min_mag1 F_min_mag1 depth_mag1
loaddatastring = '/Users/manojnair/projects/wg/MAG1Location_9data.mat';
load(loaddatastring);
load /Users/manojnair/projects/wg/HONMAY15JULY06 F_HON fday_HON

lat = lat_min_mag1(L)';
lon = lon_min_mag1(L)';
fday = fday_min_mag1(L);
F = F_min_mag1(L)';

LL = fday_HON >= min(fday) + 7/24 & fday_HON <= max(fday) + 7/24;
honf = F_HON(LL);
honfday = fday_HON(LL);
% interpolate HON data at WG sites
honfint = interp1(honfday, honf, fday + 7/24);
honfint(end) = honfint(end-1);


% interpolate HON data at WG sites
honfint = interp1(honfday, honf, fday + 7/24);
honfint(end) = honfint(end-1);
% interpolate HON data at WG sites
honfint = interp1(honfday, honf, fday + 7/24);
honfint(end) = honfint(end-1);
F_resid = ( F - nanmean(F)) - (honfint - nanmean(honfint))';
LL = isnan(F_resid);

% scatter3(lat,lon,F_resid);
% hold on;
p = fit([lat(~LL)' lon(~LL)'],F_resid(~LL)','poly55');
% plot(p);

Z = feval(p,lat,lon);
plot(fday , F_resid,'b', 'LineWidth',1);
hold on;
plot(fday , F_resid - Z,'r', 'LineWidth',1);
set(gca,'FontSize',16);
datetick('x','keeplimits');

%% Select data when MAG1 and MAG2 were together for June 24- July 7
clear
load  /Users/manojnair/projects/wg/MAG1MAG2MAY15JULY06
load /Users/manojnair/projects/wg/HONMAY15JULY06 F_HON fday_HON
load /Users/manojnair/projects/wg/MAG1Location_10data.mat L

fday_min_mag1(~L) = [];
F_min_mag1(~L) = [];
lat_min_mag1(~L) = [];
lon_min_mag1(~L) = [];
depth_mag1(~L) = [];


[fday_12,ia,ib] = intersect(fday_min_mag1, fday_min_mag2);


fmag1 = F_min_mag1(ia);
fmag2 = F_min_mag2(ib);
lat1 = lat_min_mag1(ia);
lon1 = lon_min_mag1(ia);
lat2 = lat_min_mag2(ib);
lon2 = lon_min_mag2(ib);
depth1 = depth_mag1(ia);
depth2 = depth_mag2(ib);

%limit the analyisis to days between Jule 24 to July 7

L = fday_12 > datenum(2013,6,23);

fmag1(~L) = [];
fmag2(~L) = [];
lat1(~L) = [];
lon1(~L) = [];
lat2(~L) = [];
lon2(~L) = [];
depth1(~L) = [];
depth2(~L) = [];
fday_12(~L) = [];
%remove a surevy outisde the plane

L = lon2 < -159.8518;

fmag1(~L) = [];
fmag2(~L) = [];
lat1(~L) = [];
lon1(~L) = [];
lat2(~L) = [];
lon2(~L) = [];
depth1(~L) = [];
depth2(~L) = [];
fday_12(~L) = [];

LL = fday_HON >= min(fday_12) + 7/24 & fday_HON <= max(fday_12) + 7/24;
honf = F_HON(LL);
honfday = fday_HON(LL);
% interpolate HON data at WG sites
honfint = interp1(honfday, honf, fday_12 + 7/24);
honfint(end) = honfint(end-1);

% Correct the crustal field at MAG1
p_c = fit([lat1 lon1], (fmag1),'poly11');
Z_c = feval(p_c,lat1,lon1);


% Correct MAG2 data using MAG2 data and fit the plane
p_m = fit([lat2 lon2], (fmag2-( fmag1 - Z_c) ),'poly55');
Z_m = feval(p_m,lat2,lon2);

% calculate rms error
rms1 = rms((fmag2-( fmag1 - Z_c)) - Z_m);

%plot plane & scatter
% scatter3(lat2(1:5:end),lon2(1:5:end),fmag2(1:5:end)- ( fmag1(1:5:end) - Z_c(1:5:end)) );
% hold;
% plot(p_m);



%plot time series of residual
%figure(2);
plot(fday_12, (fmag2-( fmag1 - Z_c)) -Z_m,'r');
set(gca,'FontSize',16);
datetick('x','keeplimits');
hold on;

L = isnan(honfint);

% Correct MAG2 data using HON data and fit the plane
p_h=fit([lat2(~L) lon2(~L)], ((fmag2(~L)-nanmean(fmag2(~L))) - (honfint(~L) - nanmean(honfint(~L)))),'poly55')
Z_h = feval(p_h,lat2,lon2);

% calculate rms error
rms2 = rms( ((fmag2(~L)-nanmean(fmag2(~L))) - (honfint(~L) - nanmean(honfint(~L)))) - Z_h(~L));


%plot time series of residual
plot(fday_12(~L), ((fmag2(~L)-nanmean(fmag2(~L))) - (honfint(~L) - nanmean(honfint(~L)))) - Z_h(~L), 'g')

% Consider only fmag2 and construct 




